Deep Levels and Electron Paramagnetic Resonance Parameters of Substitutional Nitrogen in Silicon from First Principles

Nitrogen is commonly implanted in silicon to suppress the diffusion of self-interstitials and the formation of voids through the creation of nitrogen–vacancy complexes and nitrogen–nitrogen pairs. Yet, identifying a specific N-related defect via spectroscopic means has proven to be non-trivial. Activation energies obtained from deep-level transient spectroscopy are often assigned to a subset of possible defects that include non-equivalent atomic structures, such as the substitutional nitrogen and the nitrogen–vacancy complex. Paramagnetic N-related defects were the object of several electron paramagnetic spectroscopy investigations which assigned the so-called SL5 signal to the presence of substitutional nitrogen (NSi). Nevertheless, its behaviour at finite temperatures has been imprecisely linked to the metastability of the NSi center. In this work, we build upon the robust identification of the SL5 signature and we establish a theoretical picture of the substitutional nitrogen. Through an understanding of its symmetry-breaking mechanism, we provide a model of its fundamental physical properties (e.g., its energy landscape) based on ab initio calculations. Moreover by including more refined density functional theory-based approaches, we calculate EPR parameters (↔g and ↔A tensors), elucidating the debate on the metastability of NSi. Finally, by computing thermodynamic charge transition levels within the GW method, we present reference values for the donor and acceptor levels of NSi.


Introduction
Nitrogen is introduced in silicon to enhance its mechanical strength [1,2], augment the precipitation of oxygen [3], suppress the diffusion of self-interstitials [3], anchor dislocations [4] and restrain the formation of voids [5]. Even if for most of these applications the precise role of nitrogen is unclear, the suppression of the formation of voids has been directly correlated to the creation of nitrogen-vacancy complexes [6]. Indeed, several nitrogen-related defects (nitrogen-nitrogen pairs, nitrogen-vacancy complexes, etc.) have been associated with the electrical and/or optical signature of N-doped silicon [7][8][9][10][11][12][13][14][15][16]. However, it has been challenging to assign a spectroscopic signal to a specific defect. In the case of the deep-level transient spectroscopy (DLTS), one set of deep levels is often assigned to the overall presence of nitrogen or to a subset of N-related defects. This ambiguous assignment can even include defects that do not share a similar atomic structure, such as susbtitutional nitrogen or the nitrogen-vacancy complex, for instance [17]. On the other hand, from electron paramagnetic resonance (EPR) spectroscopy several signals have been associated with the presence of nitrogen [18,19]. Among these EPR signals, only the so-called SL5 center has been unequivocally associated with the substitutional nitrogen (N Si ) [18][19][20][21]. In this work, we build upon this solid correlation and we establish a theoretical reference for the spectroscopic signals of substitutional nitrogen in silicon.
From the first EPR measurements [18,19], it was inferred that contrary to other group V substitutional impurities in silicon (e.g., phosphorous and arsenic), nitrogen does not remain in its tetra-coordinated on-center (T d ) site ( Figure 1). Instead, it resides in a tricoordinated off-center configuration along one of the 111 crystal directions, with a lower C 3v symmetry. This symmetry-breaking (SB) mechanism was first associated with a Jahn-Teller effect [18,19] that predicts the instability of high symmetry configurations with partially occupied degenerate states [22]. This hypothesis was rejected by early electronic structure calculations [23,24] that located the unpaired electron in a single-particle nondegenerate a 1 state of the T d configuration. The triply degenerate states t 2 were found within the conduction band and the SB mechanism was referred to as a pseudo Jahn-Teller effect [24].  The proposed symmetry-breaking (SB) mechanism T d → C 3v for this point-defect is depicted for nitrogen (in blue) and its first silicon neighbours. In the lower C 3v symmetry, the Si atom positioned along the axes of the distortion or the 111 direction is explicitly labelled. The suggested PES corresponding to this SB is also illustrated. Two equivalent off-center minima with C 3v symmetry are represented in green for the crystal directions [111] and [111]. A metastable on-center minimum (T d symmetry) is depicted in purple.
The energy landscape associated with this pseudo Jahn-Teller mechanism was investigated from electronic structure calculations by displacing the N atom along one of the equivalent 111 crystal directions [24]. The resulting potential energy surface (PES) presented a small energy barrier of ∼0.05 eV between the on-center and off-center configurations, indicating the possible metastability of the defect (see PES in Figure 1). The authors of this study, however, express caution regarding the accuracy of this value due to the simplicity of their model. Through early EPR stress measurements [18,19], the reorientation between equivalent off-center configurations at finite temperatures (T > 35 K) was observed. As a consequence of this thermally activated mechanism, the center was interpretedto be in an average tetrahedral configuration. Later EPR measurements [21,25,26] did, however, associate the existence of an on-center configuration with an increase of the hyperfine splitting when the temperature rises. This evidence was then correlated to the model of the PES in Figure 1 [21,25,26], postulating that the finite temperature dynamics of N Si are characterised by the existence of a metastable on-center geometry and the reorientation between equivalent off-center configurations.
Even if the symmetry-breaking mechanism depicted in Figure 1 was undoubtedly associated with the existence of deep levels in the silicon band gap, early optical measurements were incapable of providing a clear picture of their nature and relative position [18,19]. Furthermore, a majority of DLTS activation energies have been attributed to the overall presence of nitrogen in the sample rather than to a specific center [7][8][9][10][11][12][13]. Moreover, pre-vious studies have referred to both the substitutional impurity and the nitrogen-vacancy complex as the NV center, leading to confusion in assigning the observed levels [2,6,17,27]. Among these investigations, two independent set of deep levels were tentatively assigned to substitutional nitrogen: the first set of DLTS activation energies were attributed to a single donor and single acceptor levels at 0.31 eV and 0.08 eV from conduction band, respectively [28], whereas the second DLTS spectra were associated with single and double acceptor levels at 0.64 eV and 0.34 eV from conduction band [17].
On the computational side, recent works have performed state-of-the-art electronic structure calculations on various N-related defects in silicon [29][30][31]. While these studies have addressed the structural properties of substitutional nitrogen [29][30][31], as well as its single-particle band structure [30,31], no attempts were made to explore the possible metastability of the center or estimate its deep levels with a sufficient level of theory. In this work, we perform a thorough investigation of the SB mechanism in substitutional nitrogen from ab initio calculations, enabling us to establish its fundamental physical properties. Moreover, by including more refined DFT-based methods (i.e., linear response treatments and many-body corrections), we are able to establish a direct comparison between the atomistic model of substitutional nitrogen and the spectroscopic signals assigned to it, i.e., EPR parameters and DLTS levels.

Materials and Methods
Ab initio calculations were performed with the Quantum-ESPRESSO package for electronic structure calculations [32,33]. We used the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional [34], together with norm-conserving Trouiller-Martins pseudopotentials with GIPAW reconstruction [35,36]. We also employed a plane wave basis set with a kinetic energy cutoff of 84 Ry. A substitutional nitrogen atom is embedded in a silicon cell containing 216 atoms, which is sampled at Γ. For the neutral defect, spinpolarisation is considered. At charge states +1, −1 and −2, a compensating background charge is included.
Defect geometries at different charge states have been optimized, maintaining a constant lattice parameter. The value of this parameter was set to the estimated one for pristine silicon, a PBE = 5.47 Å. The convergence threshold on forces for the ionic minimisations was set to 0.026 eV/Å. The energy surface is explored by the climbing-NEB method [32,37] as implemented in QE. The convergence threshold of the force orthogonal to the path was fixed to 0.05 eV/Å. The hyperfine, ← → A and ← → g coupling tensors of the paramagnetic defect were computed using the QE-GIPAW module of the QE package [33,38]. This is achieved through a linear response-based approach, as described in [38,39]. To ensure convergence of these properties [40], the 216-atom supercell was integrated using a 3 × 3 × 3 Monkhorst-Pack grid [41,42].
Many-body corrections were computed within the GW method (G 0 W 0 as implemented in the ABINIT code [43,44]). Calculation parameters were kept the same as for QE, whereas the ONCVPSP pseudopotentials were employed instead [45]. Both sets of pseudopotentials showcased negligible deviation with respect to total energies. We used the Godby-Needs plasmon-pole model and a cutoff energy of 82 eV to describe the dielectric matrix. In order to assure convergence of the GW exchange-correlation self-energy, 2600 bands were employed, corresponding to an energy cutoff of 30 eV.
The thermodynamic Charge Transition Levels (CTLs) are defined as the Fermi energy at which two defect charge states have the same formation energy [46]. In simpler terms, a CTL represents the energy necessary for a transition to take place between two different charge states of a defect. In the particular case of substitutional nitrogen, such a transition can be expressed as, where N q Si and N q Si represent the defect at charge states q and q . If we consider the formation energy E q f of the charge state N q Si as the reference, the formation energy of the charge state N q Si can be expressed as, In this equation, ∆(q/q ) represents the energy exchanged during the transition between charge states q and q . This difference accounts for the energy required to add or remove an electron from N q Si , along with the energy released due to atomic relaxation. The chemical potential of electrons, µ e , is often written as the Fermi energy referred to the valence band maximum (VBM), , where CBM is the conduction band minimum. The CTL is then simply, From Equation (3), we evaluated the CTLs using two different approaches: the standard PBE approximation and the so-called DFT and GW combined approach [47]. Within the first method, ∆(q/q ) is simply evaluated as the difference in DFT total energies between the two relaxed defect structures at charge states q and q , where R q denotes the ground state geometry at charge state q. The ionization potential of the VBM, VBM , is approximated by the corresponding Kohn-Sham eigenvalue in pristine silicon. In order to circumvent the self-interaction contribution in standard PBE calculations, the second approach relies on the GW method to determine accurate ionization potentials (IPs) and electronic affinities (EAs) [48,49]. The difference ∆(q/q ) is thus determined as an IP/EA, evaluated as a quasiparticle (QP) energy at charge state q, minus/plus a relaxation energy. This relaxation energy is evaluated as a DFT total energy difference between the relaxed structures at q and q , evaluated for the final q charge state. The ∆(q/q ) difference is thus evaluated as, Finally, VBM in Equation (3) is evaluated as a QP energy of pristine Si. The formation energy of charged defects (Equation (2)), is often corrected to mitigate the spurious electrostatic interaction between defect replicas [50]. In the present work, the correction is evaluated by aligning the Density Of States (DOS) of the charged and neutral defect states. In the case of the positive charge state, this shift is equal to −0.03 eV, whereas for both singly and doubly negative charge states it is equal to 0.01 eV.

The Symmetry Breaking in N Si
Our first-principle calculations confirmed that the ground state geometry of N 0 Si displays an off-center configuration with C 3v symmetry, rather than the on-center configuration or the tetrahedral (T d ) symmetry. In the on-center configuration, four equivalent N-Si interatomic distances of 2.05 Å are found between the N atom and its nearest neighbours (or silicon atoms 1, 2, 3 and 4 in Figure 1). These Si atoms form a tetrahedron with an edge length of 3.34 Å. With the atomic distortion, the N atom moves away from one of its nearest neighbours (silicon atom 4 in Figure 1), elongating the N-Si distance to 3.15 Å (Figure 2a). The remaining three equivalent N-Si distances are shortened to 1.86 Å, similar to the N-Si bond length (1.73 Å) in silicon nitride [51]. The displacement of the N atom along the C 3 symmetry axis also distorts the tetrahedron formed by its nearest neighbours. The Si atom along the axis of distortion (denoted from now on the Si 111 atom) is 3.87 Å away from the remaining three Si atoms, whose interatomic distances are 3.20 Å. Our analysis of the symmetry breaking indicates that the total displacement of the nitrogen atom away from its substitutional site is 0.57 Å. In its off-center configuration, the N atom lies relatively close (0.25 Å) to the interatomic plane defined by its nearest neighbours (or silicon atoms 1, 2 and 3 in Figure 1). Both these distances are in agreement with previously reported values [30,31]. The DFT total energy is plotted with respect to the displacement of the N atom from its on-center position to its off-center configuration. The energy difference between the two minima is shown, as well as the energy barrier from the off-center to the on-center geometries.
To gain a deeper understanding of the symmetry breaking (SB) in substitutional nitrogen, we examined the one-electron defect states of both the on-center and off-center configurations. We approximated the one-electron wavefunctions using Kohn-Sham eigenfunctions and their corresponding energy levels using quasiparticle (Kohn-Sham) eigenvalues computed at Γ. Figure 2b shows the one-electron energies corresponding to the on-center and off-center geometries, along with the real space projection of the wavefunctions of the singlet defect states. In the on-center configuration, we identified a partially occupied a 1 state at 0.43 (0.20) eV from the valence band maximum (VBM). As depicted in Figure 2b, the one-electron wavefunction is centered in the N atom and is equally distributed among the four neighboring Si atoms. The unoccupied triplet states t 2 , predicted in the theoretical model [24] (Figure 1), are found within the conduction band, at 0.09 (0.08) eV from the conduction band minimum (CBM). After the symmetry breaking, the partially occupied a 1 state lowers its energy, whereas the t 2 state splits into a second singlet a 1 and a doublet e (in agreement with symmetry considerations [52] and the previous model [24]). The partially occupied singlet state is located within the band gap, at 0.07 (0.02) eV from the VBM, and is primarily localized on the Si 111 atom and to a lesser extent on the nitrogen atom. The second singlet state is found within the conduction band, at 0.03 (0.03) eV from the CBM. The doublet state hybridizes with the conduction states and is found at 0.11 (0.11) eV from the CBM. Notice that both the KS wavefunctions and the QP eigenvalues are computed at Γ and represented within the folded bulk bands. When k-points are added to the calculation, we, however, observe a large dispersion of the defect states that is more prominent for the a 1 state of the on-center configuration (for both spin channels) and the empty a 1 state of the off-center configuration (minority spin channel). The band structures for both defect configurations are shown in the Supplementary Materials. Based on the energy diagrams shown in Figure 2b and, in particular, the absence of partially occupied degenerate states, we confirm that the symmetry breaking T d → C 3v is not a Jahn-Teller effect [22]. Thus, the PES of substitutional nitrogen does not exhibit the standard conical intersection, where the highly symmetric configuration is located at the top of the so called Mexican hat-shaped PES. In this case (or in the case of a pseudo Jahn-Teller effect [24]), the shape of the PES cannot be induced from symmetry considerations and so the nature of the critical points (T d , C 3v ) of the PES cannot be predicted without explicit electronic structure calculations. To begin exploring the PES of substitutional nitrogen in silicon, we examined the Minimum Energy Path (MEP) between the on-center and off-center configurations. In Figure 2c, the total DFT energy is plotted along the MEP, whose coordination reaction is taken as the displacement of the N atom away from its ideal on-center position. As depicted in Figure 2c, the MEP between these geometries corresponds to a linear interpolation between both structures. It should be noted that the energy difference between the two configurations is not negligible, with a value of 81 meV. This value is in decent agreement with the previous estimation of 50 meV [30] and 80 meV [29]. Additionally, we observe an energy barrier of 129 meV between the on-center and off-center configurations. Our results therefore indicate that both the tetrahedral (T d ) and the trigonal (C 3v ) configurations are metastable.
We further explored the PES of substitutional nitrogen by investigating the transition path between equivalent off-center configurations. Since the symmetry breaking preserves any of the original C 3 symmetry axes, the nitrogen atom can move along any of the equivalent 111 directions away from its first silicon neighbors ( Figure 1). As a result, the PES is characterised by four equivalent global minima as off-center configurations and the defect can reorient at finite temperatures. As a first approximation to the MEP between equivalent configurations, we tested a linearly interpolated path along which the N atom avoids its on-center site. However, the transition path relaxes to pass through the tetrahedral (T d ) symmetry, resulting in the MEP in Figure 2c. The N atom is therefore forced to transit through its on-center configuration to jump between equivalent off-center configurations. It is important to note that the activation energy for this reorientation mechanism cannot be simply estimated as the energy barrier between the on-center and offcenter configurations, since the T d symmetry is a metastable minimum. Namely, the defect could remain trapped in its on-center configuration, delaying the transition to another global minimum. Before comparing our shape of the PES to experimental observations, we address the potential influence of size effects on the MEP. For this purpose, we performed relaxation calculations for both on-center and off-center defect structures within a 512-atom cell. The resulting interatomic distances differ by no more than 0.03 Å from the previously described bond lengths. The energy difference between these structures is 81 meV and the energy barrier is 129 meV, indicating that the MEP obtained with the 216-atom cell converged within 10 meV.

Ab Initio EPR Parameters
Once the ground state configuration of the neutral substitutional nitrogen is well described, we calculate its corresponding EPR parameters to provide an explicit comparison with experimental data. In particular, in Table 1, we show, for the off-center configuration, the principal values of the ← → g tensor and of the hyperfine tensors together with their isovalues (g iso and A iso ). The hyperfine tensors refer to the silicon atom along the atomic distortion (or the Si 111 atom) and to the nitrogen impurity. The substantial disparity between A i principal values corresponding to the Si 111 atom and the N atom reflects the notable localization of the unpaired electron on the Si atom (see wavefunction of the a 1 state in Figure 2b).  [18][19][20]. Calculated hyperfine ← → A tensors are also displayed for the on-center geometry in the 216-atom cell (tw-on) and in the NSi 4 H 12 cluster (tw-cluster). These values are compared to EPR data [21,26]. Also, the ← → g tensor as calculated for the NSi 4 H 12 cluster is given. Concerning the ← → g tensor, we account for the axial symmetry of the defect and we designate the principal value g 1 as the one corresponding to the principal axis parallel to the 111 direction (or g ). Consequently, the g 2 and g 3 values correspond to any two axes perpendicular to the 111 direction (or g ⊥ ). From our values in Table 1, it is apparent that the axial symmetry is not fully preserved, as indicated by the inequality of g 2 and g 3 . This small deviation is within the expected inaccuracy estimates of the GIPAW method [53,54] and cannot be attributed to a poor description of the defect geometry, since the estimated interatomic distances are fully compliant with the C 3v symmetry. However, it can be assigned to the artificial dispersion of the defect states when several k-points are included in the sampling of the Brillouin zone (see band structures in the Supplementary Material).
In Table 1, we compare the EPR parameters calculated for the off-center configuration to the reported EPR data of the SL5 center [18][19][20]. The calculated principal values of the hyperfine structure tensor ← → A are within 12% of the experimental ones for the silicon atom along the 111 distortion and within 20% for the nitrogen impurity [18][19][20]. Moreover, both the principal values and the isovalue of the ← → g tensor lie within ∼1000 ppm of those measured via Brower [18,19] and with even better agreement (∼700 ppm) with the most recent measurements [20].
Earlier theoretical calculations of the hyperfine tensors, for the off-center configuration, were derived using the Hartree-Fock approximation for a relatively small atomic cluster NSi 4 H 12 [55]. However, the calculated isovalues exhibited significant discrepancies when compared to experimental EPR data [18][19][20]. Specifically, the calculated value for A iso ( 29 Si) was overestimated, while the one for A iso (N) was underestimated, with deviations of nearly a factor of 2 in both cases [55]. By contrast, our ab initio EPR parameters (both ← → g and ← → A tensors) display excellent agreement with the experimental data and hence provide a firmly grounded theoretical assignment of the SL5 EPR center to the distorted off-center configuration of the substitutional nitrogen.
Besides investigating the magnetic properties of a frozen N atom, the dynamics of the SL5 center were also investigated from EPR stress measurements [18,19]. In particular, the broadening of the spectral lines at temperatures above 35 K was attributed to the reorientation of the defect between equivalent off-center configurations. Based on this analysis, the activation energy for the N atom to jump from one 111 distortion to another was inferred to be 0.11 ± 0.02 eV [19]. From our exploration of the PES, we calculate an energy barrier of at least 0.129 eV between equivalent off-center distortions (Figure 2c). This evaluation is in rather good agreement with the experimental estimate of 0.11 eV. Notice, however, that given the existence of a metastable on-center minimum in our reorientation path, the effective activation energy for the reorientation process is expected to be slightly higher. Before even estimating this effective barrier, in the following, we reassess the existence of a metastable on-center configuration.
From EPR measurements [18][19][20][21]25,26], an argument in favor of a metastable oncenter configuration was put forward to explain the observed increase in the hyperfine splitting when augmenting the temperature [21,26]. In order to elucidate the existence of this metastable minimum, we estimated its EPR parameters (Table 1) and compared these values to the ones previously assigned to this tetrahedral configuration [21,26]. In this highly symmetric configuration, the ← → g tensor and the ← → A (N) hyperfine tensor are expected to be isotropic. In fact, as reported in Table 1, we found an isotropic hyperfine tensor with a considerably large isovalue (165.7 MHz), while the isovalue on the silicon (|A iso ( 29 Si)| = 40.1 MHz) is about the same as found for the nitrogen impurity in the off-center configuration. Because of numerical instabilities in the calculations of the ← → g tensor for the on-center configuration, we performed an EPR demonstration calculation on a hydrogen-passivated silicon cluster; specifically, the NSi 4 H 12 cluster. This cluster features four N-Si bonds with a length of 1.89 Å and Si-Si distances of 3.08 Å. As shown in Table 1, both the calculated ← → g and ← → A (N) tensors are isotropic and the value of the hyperfine coupling constant (226.6 MHz) is consistent with the one (165.7 MHz) calculated for the on-center N Si in the silicon supercell.
The calculated hyperfine tensors ← → A (N) of the on-center configuration are in disagreement with the experimental room temperature hyperfine splitting value given in [21]. Yet, we note that in [26] the temperature dependence of the hyperfine splittings is fitted by means of Boltzmann weights. From the fit, a remarkably large value of 376 MHz was assigned to the N in the on-center configuration. From our model, the defect state is transferred from N to Si when moving from on-center to off-center configuration (see Figure 2c). This implies that the value of A iso (N) drastically decreases with this transfer; thus, it would seem unlikely that this value only differs by about 5 MHz between the two configurations [19,21,28]. A remarkably larger difference can be expected, as suggested by Murakami's fit [26] and by our results presented in Table 1. In fact, A iso should be significantly larger for the atom around which the paramagnetic state is located. In our calculations (Table 1), we see, for instance, a A iso ( 14 N) and a A iso ( 29 Si) of 226.6 MHz and −28.3 MHz, respectively, for the on-center configuration that become 32.2 MHz and −264.1 MHz, respectively, in the off-center configuration. From our computational model, we remark that even if a metastable on-center minimum is found in the exploration of the PES, the ratio between the off-center and on-center occupation (from Boltzmann weights) is about 120 at room temperature. Hence, its implication in the thermally activated dynamics is minimal and its contribution to the averaged spectroscopic signals is minor, leading just to a small increase (∼5 MHz) in the hyperfine splittings at room temperature. Yet, the on-center configuration can be tracked down from the temperature-induced increase in the hyperfine splittings [26]. The latter, however, might also involve an effect due to coupling with the phonons of the off-center potential well.
We note that, contrary to the interpretation suggested by Ref [21], the on-center configuration should not be mistaken with the room-or high-temperature averaged configuration of N Si . In fact, the averaged spectroscopic signals depend on the occupation probability of both minima (the four degenerate off-center and the on-center).

Charge Transition Levels of N Si
Once a computational model for neutral substitutional nitrogen has been proposed, we investigate its relative stability compared to other charge states. In particular, we investigate the charge states +1, −1 and −2 and we determine their respective CTLs. At charge state +1, there is no trapped electron within the silicon band gap and the defect remains in an on-center configuration, with a N-Si bond length of 2.05 Å. At charge state −1, the offcenter configuration is preferred and so the a 1 defect state is fully occupied. Along the distortion, the N-Si 111 bond is 3.07 Å, whereas the remaining three N-Si distances are equal to 1.85 Å. If yet another electron is added to the defect, the system remains in an off-center configuration, with a characteristic N-Si 111 distance of 1.85 Å, compared to the three N-Si bond lengths of 3.04 Å. This charge state −2 is paramagnetic, with the defect state a 1 in the conduction band partially occupied.
We present our calculated thermodynamic CTLs, obtained within the PBE+GW (PBE) approach, for the single donor and the single acceptor levels of substitutional nitrogen in Table 2. The donor level is found relatively close to the valence band maximum, with calculated values of VBM + 0.19(0.08) eV. In contrast, the acceptor level is located deep in the silicon band gap, positioned at CBM − 0.55(0.32) eV from the conduction band minimum. Notice that both values correspond to the mid-region of the silicon band gap, considering that the band gap is estimated to be 1.02 (0.63) eV. Finally, the double acceptor level is found within the conduction band, at 0.06 (0.06) eV from the conduction band minimum. The charge state −2 is therefore thermodynamically unstable.
The existence of a single donor and acceptor levels for N Si was previously reported via an ab initio calculation [29]. The proposed values were VBM + 0.5 eV and CBM − 0.4 eV, respectively. Notice, however, that they reported an electrical gap of 0.5 eV and that a corrective shift of the estimated levels was made. Contrary to our approaches (Table 2), their corrected levels are both deep in the band gap, at ∼0.25 eV from each other. Given the well-established accuracy of the DFT + GW approach estimating single donor and acceptor levels, we emphasize the existence of a donor level close to the valence band maximum.
In Table 2, we present a comparison between our CTLs and the DLTS activation energies tentatively assigned to the off-center substitutional nitrogen [17,28]. We observe a significant disparity between the reported values in [28] and our estimated levels. Given the accuracy of the DFT + GW method [47], we argue that the assignment of the levels CBM − 0.31 eV and CBM − 0.08 to substitutional nitrogen was incorrect. In contrast, we notice the proximity of our estimated acceptor level with the one proposed in a recent DLTS investigation [17]. Nonetheless, we acknowledge the discrepancy between their identified double acceptor level and our computed value. From this disparity, we conclude that either the assignment of both acceptor levels to substitutional nitrogen was incorrect or that these levels are the signature of another N-related defect (e.g., the nitrogen-vacancy complex). Table 2. Deep levels located within the bandgap associated with substitutional nitrogen in silicon. Thermodynamic CTLs were calculated within a standard PBE approximation (tw-PBE) and within the combined DFT and GW approach (tw-PBE+GW). Previously estimated levels (pw) are also included [29]. DLTS levels were tentatively assigned to susbstitutional nitrogen [28] and activation energies associated with the off-center nitrogen, however, also addressed as NV center, in [17]. All levels are given from the conduction band minimum in eV.

Conclusions
We performed a systematic investigation of the symmetry-breaking mechanism in substitutional nitrogen from first-principle calculations. We characterized fundamental properties of both the on-center and off-center configurations of N Si , including the relative position of its single-particle defect states within the band gap. Moreover, we proposed a model for its potential energy surface by performing an extensive search of the minimum energy path between equivalent off-center configurations. Such a transition path forces the nitrogen center to pass through the on-center configuration that is estimated to be a shallow metastable minimum.
In order to provide a direct comparison with spectroscopic signals, we computed the EPR parameters and the thermodynamic CTLs of N Si . Our EPR values for the offcenter configuration are in excellent agreement with the SL5 EPR signature. Moreover, the energy barrier we calculate for the reorientation mechanism between equivalent off-center configurations is in fair agreement with EPR stress measurements, validating our model for its potential energy surface. Furthermore, our calculations further support a large hyperfine splitting value for the N Si in the on-center configuration as estimated in [26]. Finally, we proposed a set of deep levels for N Si , corresponding to a single donor and acceptor levels. Our values allowed us to identify incorrect assignments between DLTS activation energies and the electrical activity of N Si and they can guide its identification in future works.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/nano13142123/s1. Supplementary materials include data on band structure of the defective silicon cells ( Figure S1).